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Ground-state solutions in a dilute gas interacting via con- 
tact and magnetic dipole-dipole forces are investigated. To 
the best of our knowledge, it is the first example of studies 
of the Bose-Einstein condensation in a system with realis- 
tic long-range interactions. We find that for the magnetic 
moment of e.g. chromium (6/zb) and a typical value of the 
scattering length all solutions are stable and only differ in 
size from condensates without long-range interactions. By 
lowering the value of the scattering length we find a region 
of unstable solutions. In the neighborhood of this region the 
ground state wavefunctions show internal structures not seen 
before in condensates. Finally, we find an analytic estimate 
for the characteristic length appearing in these solutions. 



Since the advent of Bose-Einstein condensation in di- 
lute gases of alkalies [0 and hydrogen Q it has become 
apparent that the interactions between the condensed 
atoms govern most of the observed phenomena. In all 
the experiments so far the interaction can be described 
by a contact potential which is characterized by the scalar 
quantity a being the s-wave scattering length. Static 
properties like the condensate's ground-state density pro- 
file, its instability in the case of negative a and the mean- 
field shift in spectroscopic measurements as well as dy- 
namic properties like collective excitations and propaga- 
tion of sound have been investigated ||] . Similarly, non- 
linear atom optics experiments e.g. four wave mixing Q] 
are only possible due to the large nonlinearity mediated 
by the atom-atom interactions. The T=0K situation in 
almost all experiments can be described very well by the 
Gross-Pitaevskii equation (s| . 

Any reasonably strong dipolc dipolc interaction would 
largely enrich the variety of phenomena to be observed 
in dilute gases due to their long range and vectorial char- 
acter. However, for all of the condensed atomic species 
the magnetic moment /z was roughly 1 /is (Bohr mag- 
neton) and the respective magnetic dipole-dipole inter- 
action was negligible compared to the contact potential. 
It has been proposed to induce a strong electric dipole- 
dipole interaction in alkalies by the application of strong 
DC electric fields §. 

Recently it has become possible to trap atoms with 
higher magnetic moments at high densities. Examples 
are europium (/i = 7/is) 0, which has been trapped 
magnetically by a buffer-gas loading technique, and 
chromium (/i = 6(J-b), which has been loaded into a mag- 



netic trap by a buffer-gas technique Q and by laser cool- 
ing |)| . For these species the scattering lengths are not 
known to date, but, assuming a normal non-resonant be- 
havior, the crossections of the scalar contact potential 
and the magnetic dipole-dipole interaction are of com- 
parable size. 

As the dipole-dipole interaction is attractive parallel 
to a common polarization axis the immediate question 
arises: can a stable condensate be formed under the in- 
fluence of a dipole-dipole interaction? What is its effect 
on anisotropic clouds? What do the ground-state wave- 
functions look like? In this Letter we address these ques- 
tions and choose the magnetic dipole-dipole interaction 
in a dilute cloud of chromium atoms as an example for 
numerical calculations. The results however apply for all 
static dipole-dipole interactions. 

In general, the Gross-Pitaevskii equation with a long- 
range interaction term and for a cylindrical harmonic 
trap has the following form: 
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where ^ is the mean-field condensate wavefunction, a - 
the s-wave scattering length, N - number of atoms and 
m - mass of the atom. The reference frequency of the 
trap wo is chosen in the xy plane and the anisotropy of 
the trap is defined by the 7 factor. In our case, the long- 
range potential takes the form characteristic of magnetic 
dipole-dipole interactions: 

= Mo l?i{r) ■ m(r') - 3/gXg) ■ u py(P) ■ u 

where u = S~^Ti ana - Mo is magnetic permeability of the 
vacuum. We will assume that all the magnetic moments 
are in the same direction (/Ti = /I2 ) , which will be referred 
to as the polarization direction. Note that, depending on 
a configuration of dipoles, magnetic potential can be re- 
pulsive as well as attractive. This fact is a source of a 
variety of phenomena which do not appear in conden- 
sates with contact interactions only. Another peculiar 
feature of this long-range (integral) term is that for a uni- 
form density distribution it vanishes if integrated within 
a sphere, which means that in that case a dipole placed 
in the middle of the sphere would feel no magnetic force. 
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Equation (Q) is now an integro-differential equation. 
The integral term can be simplified, because a part of it 
can be calculated analytically. First of all, one should 
notice that this term has a form of a convolution: 



V(f- r')\$>(r')\ 2 d 3 r' = V{r) * |#(f)| 2 



The transform of the potential reads: 

-cnrt^w 2/i q 2 s r coa(bq) sin(6g) 

J-{V[r)) = nofi (1 — 3 cos a) 



(bq) 2 (bq) 3 



(3) 
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where b is a distance below which atoms overlap (i.e. the 
radius of the atoms which is of the order of a few Bohr 
radii) and a is the angle between the Fourier variable 
q and ft. A value of b is small in comparison with a 
length scale in the system (oscillator unit) and we will 
not consider large values of q therefore it is sufficient to 
use the limit: 

lim T{V{r)) = -\^ 2 (1 - 3cos 2 a). (5) 

b^O 3 

In order to evaluate the Fourier transform of eq. (3) 
JFdV'l 2 ) is calculated numerically (FFT) and multiplied 
by J-(V(r))which depends solely on a. 

In order to obtain the ground-state solution for the 
condensate one has to solve the Gross-Pitaevskii equation 
in the form of Eq. (|l|) . To do it we employ a variant of the 
split-operator method called the imaginary-time propa- 
gation method, being now a common routine of solving 
the Gross-Pitaevskii equation. For each time-step one 
needs to compute four FFT's: two for the long-range 
term calculation and two for the evolution. 

Now we will list parameters used in the calculations. 
The reference trap frequency ujq is equal to 2n 150 Hz 
which is an accessible and typical value in present exper- 
iments. This corresponds to a characteristic length unit 

in the system (oscillator unit) d = J m ^ uq ~ 1.14 /im. 

The s-wave scattering length for chromium is unknown 
- instead a value for sodium is assumed tentatively 
a = o-Na = 2.75 nm. 

One should keep in mind that all the presented so- 
lutions are scalable. For the Gross-Pitaevskii equation 
with as well as without long-range forces the solutions 
stay the same - in scaled coordinates - as long as the 
product N^/ujq is kept constant. 

The ratio of the contact and long-range terms in gen- 
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eral is given by ^— — , m being the atom's mass, and solu- 
tions stay identical if it is kept constant. For the chosen 
example of chromium we will in the second part of the 
Letter vary the scattering length to change that ratio. 

The main conclusion is that, for this set of param- 
eters (in particular the scattering length value) all the 
investigated solutions are stable, even in very asymmet- 
ric traps. Moreover, the wavefunctions possess familiar 
shapes: Gaussian-like for a small number of atoms and 
parabolic- like for bigger condensates. 



At first, we investigated two extreme cases in cylindri- 
cal traps. Following simple intuition about an interac- 
tion of two dipoles, we designed these situations so that 
in the first one there would be a majority of repulsive 
forces, whereas in the other one an attractive component 
would be dominant. 

In the first situation, the trap is flattened in a plane 
perpendicular to the polarization axis (xy plane) with 
the asymmetry factor 7 = ^ — 10 (referred to as disc- 
shaped). For this, and only this case, we lower the refer- 
ence frequency to wq = 2n 15 Hz to make tu z = IOwq 
experimentally accessible (the corresponding oscillator 
unit is d =~ 3.6 fj,m). In the inset of Fig|l] half- width 
of the ground state vs. number of atoms is plotted for 
the disc-shaped trap. It is compared to a corresponding 
half-width for a Thomas-Fermi solution for a condensate 
with contact interactions only (referred to as the stan- 
dard Thomas-Fermi limit). Note that the Thomas-Fermi 
approximation is good only for large number of atoms, 
thus for N < 10,000 it has not been used. In the z- 
direction the exact solution is indistinguishable from the 
standard Thomas- Fermi limit, because z is a stiff direc- 
tion and so all the interparticle interaction effects are 
small with respect to the trapping forces. On the other 
hand, along the x-axis the condensate expands with re- 
spect to the standard Thomas-Fermi limit as it should be 
the case for the repulsive direction. To conclude, in the 
disc-shaped trap the net effect of interactions is repulsive 
and the condensate expands. 

The opposite extreme case is a trap stretched along 
the polarization axis (z-axis) with 7 = 0.1 (referred to as 
cigar-shaped) - see Fig|l]. For the soft direction (along 
the z-axis, being also an attractive direction), in agree- 
ment with simple intuition, one observes shrinking of the 
condensate. The differences with respect to the standard 
Thomas-Fermi limit for a condensate without long-range 
interactions diminish as the number of atoms grows. This 
is the case because in the z direction the wavefunction 
is very flat (roughly uniform) and so the integral compo- 
nent almost vanishes. Surprisingly, we observe shrinking 
of the condensate in the x (repulsive) direction as well. 
We argue that it is caused by the attraction exerted by 
very many dipoles concentrated along the soft (z) direc- 
tion which takes over the repulsive contribution along the 
x-axis. Concluding, in the cigar-shaped trap we observe 
an overall shrinking of the condensate. 

As we see, long-range magnetic interactions do not 
change qualitatively the ground-state solutions of the 
Gross-Pitaevskii equation, but they affect mainly a size 
of the condensate. This indicates that for the used set 
of parameters an influence of contact (repulsive) inter- 
actions is still dominant. However, as we do not know 
the scattering length for chromium and there are several 
methods proposed to manipulate this parameter by ap- 
plication of various external fields (e.g. see JToj] , [Tlfl), we 
decreased the assumed value and by this means enhanced 
the effect of long-range (partially attractive) forces. As 
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a result, we were able to find a region of unstable so- 
lutions. An even more striking discovery was the ob- 
servation of structured shapes of the wavefunctions ac- 
quired by the condensates near the instability thresh- 
old. By probing the parameter space (a, N) (scattering 
length, number of atoms) we found the phase diagram 
depicting localizations of stable, structured (still stable) 
and unstable solutions shown in Fig.^j. For simplicity 
we did calculations for a spherical trap with a frequency 
u>q=2tt 150 Hz, but a few cases calculated for asymmetric 
traps convinced us that the qualitative stability behav- 
ior would not change. The polarization direction is still 
z. For illustration purposes, we present two examples 
of the structured solutions. The first one (Fig.||) was 
obtained for 80,000 atoms with a/a,/v Q =0.115 which sit- 
uates this case very near the instability threshold. The 
second example (Fig.|J) is a wavefunction for 4,000,000 
atoms and a/ajv a =0.233 (again right above the instabil- 
ity threshold). We have also used geometry of the trap 
to stabilize solutions unstable in a spherically symmetric 
potential (by increasing the asymmetry factor 7) as well 
as to destabilize ones that were stable in a spherical trap 
(by decreasing 7). It is remarkable that the complex- 
ity of the ground-state wavefunction arises solely from 
atom-atom interactions and does not reflect simplicty of 
an external potential (it represents self-organization in 
the ground state). 

For all the structured solutions there exists a char- 
acteristic length defined by a distance between adja- 
cent maxima in the plane perpendicular to polariza- 
tion, which is roughly 3d. In order to explain this fea- 
ture, we performed a stability analysis for an infinite (no 
trap), homogenous case. In this case a proper solution is 
* = -i exp(-igt), V being the volume and g = N^^. 
By imposing a small harmonic perturbation of the form 

= -^[1 + u {t) cos(q ■ r)] exp(-igt) we found that the 
unstable perturbations are those for which: 

2 4JVm Cr w 2 Awha 

q < —— (1-3 cos a) , 6 

nv oh rricr 

where a is the angle between the wavevector and the po- 
larization direction. This result implies that surprisingly 
all the perturbations parallel to polarization (a! = 0) are 
stable, whereas for the direction perpendicular to polar- 
ization (a = ^ ) there is a long- wave instability. To com- 
pare this result with the observed characteristic length of 
3c? we crudely approximated the volume of our conden- 
sates by V = ^7vRtp, where Rtf is the Thomas-Fermi 
radius for a condensate without long-range interactions. 
This yields the critical wavelength: 

A 2 _ An 3 h (15a)idir 

For the two structured wavefunctions presented in this 
letter, this value corresponds to 2. 2d and 1.8d, respec- 
tively, which, for such an approximate analysis, is in a 



surprisingly good agreement. We note in passing that the 
onset of unstable perturbations (Ac — > 00) corresponds 
to ac — 0.3et7v a which, again, roughly armroximates the 
calculated instability threshold (see Fig.g). 

Motivated by experiments under development, all our 
calculations used chromium parameters. They showed 
that ground-state solutions differ in a non-trivial way 
from the usual solutions and found a region of insta- 
bility as well as structured solutions in its neighbor- 
hood. Reaching the interesting region near the instability 
boundary for chromium seems problematic with state- 
of-the-art techniques. However, in the case of europium 
(larger mass and magnetic moment) the instability 
threshold for 1,000,000 atoms is situated at 92% of the 
sodium scattering length value and still grows for larger 
numbers of atoms (e.g. 108% for 10,000,000 atoms). Our 
results may also find use in the blooming area of cold 
molecules |l2]| some of which possess large permanent 
electric dipole moments [[l3| . For a typical electric dipole 
moment value of 1 Debye a prefactor of the long-range 
term is 4 orders of magnitude greater than the corre- 
sponding term for magnetic dipole of 1 Bohr magneton. 
One can thus conclude that behavior of a polar molecular 
BEC will be overruled by dipole forces. 
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FIG. 1. Half-width of the chromium condensate vs. num- 
ber of atoms in a cigar-shaped cylindrical trap with 7 = 0.1 
and a=aMa- Filled circles indicate an exact solution for a 
section along the :r-axis as compared to the corresponding 
result in the standard Thomas- Fermi limit (empty circles). 
Empty squares represent an exact solution along the z-axis in 
comparison with the corresponding standard Thomas-Fermi 
limit (filled squares) . The inset presents analogous results in a 
disc-shaped trap (7 = 10) and for a—aNa- An exact solution 
for a section along the z-axis (filled squares) is indistinguish- 
able from the standard Thomas- Fermi limit. 

FIG. 2. Stability diagram for a spherical trap with 
u>o=2n 150 Hz. Number of atoms is plotted on the horizontal 
axis and the scattering length (expressed as a fraction of the 
corresponding value for sodium) on the vertical one. Solid 
fine is the instability threshold - below it no stable solutions 
could be obtained. The dashed line is the boundary between 
standard-shaped and structured solutions. In the shaded area 
structured stable solutions are present. Error bars indicate 
discreteness of our parameter-space probing. 
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FIG. 3. Squared modulus of the ground-state wavefunction 
in the xz plane for a spherical trap with u>o=2-k 150 Hz, 80,000 
atoms and the scattering length a/ajv a =0.115. The horizontal 
axes are in oscillator units d. 



FIG. 4. Squared modulus of the ground-state wavefunc- 
tion in the xz plane for a spherical trap with luo=2tt 150 Hz, 
4,000,000 atoms and the scattering length a/ajv a =0.233. The 
horizontal axes are in oscillator units d. 
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